clear;
T=30;
N=10000;
dt=T/N;
t=0:dt:T;

y=zeros(3,length(t));
s=10;
b=8/3;
r=28;
y(:,1)=[-1.6;-0.5;21];

for n=1:length(t)-1
    k1=[s*(y(2,n)-y(1,n));r*y(1,n)-y(2,n)-y(1,n)*y(3,n);-b*y(3,n)+y(1,n)*y(2,n)];
    k2=[s*((y(2,n)+k1(2)*dt/2)-(y(1,n)+k1(1)*dt/2));r*(y(1,n)+k1(1)*dt/2)-(y(2,n)+k1(2)*dt/2)-(y(1,n)+k1(1)*dt/2)*(y(3,n)+k1(3)*dt/2);-b*(y(3,n)+k1(3)*dt/2)+(y(1,n)+k1(1)*dt/2)*(y(2,n)+k1(2)*dt/2)];
    k3=[s*((y(2,n)+k2(2)*dt/2)-(y(1,n)+k2(1)*dt/2));r*(y(1,n)+k2(1)*dt/2)-(y(2,n)+k2(2)*dt/2)-(y(1,n)+k2(1)*dt/2)*(y(3,n)+k2(3)*dt/2);-b*(y(3,n)+k2(3)*dt/2)+(y(1,n)+k2(1)*dt/2)*(y(2,n)+k2(2)*dt/2)];
    k4=[s*((y(2,n)+k3(2)*dt)-(y(1,n)+k3(1)*dt));r*(y(1,n)+k3(1)*dt)-(y(2,n)+k3(2)*dt)-(y(1,n)+k3(1)*dt)*(y(3,n)+k3(3)*dt);-b*(y(3,n)+k3(3)*dt)+(y(1,n)+k3(1)*dt)*(y(2,n)+k3(2)*dt)];
    y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4)*dt;
    
end
subplot(1,2,1)
plot(t,y)
subplot(1,2,2)
plot3(y(1,:),y(2,:),y(3,:))
